支配方程式

---

連続の式

\begin{equation} \nabla\cdot \boldsymbol{v} = 0 \tag{1} \end{equation}

非圧縮ナビエストークス方程式

\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v}\cdot\nabla)\boldsymbol{v} = - \nabla p + \frac{1}{Re}\Delta \boldsymbol{v} \tag{2} \end{equation}

$Re$はレイノルズ数．

圧力の楕円型方程式

(2)式の発散をとって， \begin{equation} \frac{\partial D}{\partial t} + \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) = - \Delta p + \frac{1}{Re}\Delta D, \quad D \equiv \nabla \cdot \boldsymbol{v} \notag \end{equation} \begin{equation} \therefore \: \Delta p = -\frac{\partial D}{\partial t} - \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) + \frac{1}{Re}\Delta D \end{equation} $D=0$という条件を使えば， \begin{equation} \therefore \: \Delta p = - \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) \tag{3} \end{equation} というポアソン方程式を得る．しかし数値的には必ずしも$D=0$になるとは限らないから時間微分項は，数値計算するときを想定して，単純に

\begin{equation} \frac{\partial D}{\partial t} \simeq \frac{D_{n+1} - D_n}{\Delta t} \end{equation}

としておいて，次の時刻$t_{n+1}$で$D$がゼロになるように

\begin{equation} \frac{\partial D}{\partial t} \simeq -\frac{D_n}{\Delta t} \end{equation}

というふうにする．あと実際$D$は小さい数になると思われるので2階微分はほぼゼロとして消す．なお，

\begin{align} \Delta\boldsymbol{v} &=\mathrm{grad}\:\mathrm{div}\boldsymbol{v}-\mathrm{rot}\:\mathrm{rot}\boldsymbol{v} \\ \therefore \: \mathrm{div} \Delta\boldsymbol{v} &= \mathrm{div} \mathrm{grad}(\mathrm{div}\boldsymbol{v})-\mathrm{div} \mathrm{rot} \:(\mathrm{rot}\boldsymbol{v}) \\ &= \Delta (\mathrm{div}\boldsymbol{v}) = \Delta D \end{align}

を用いた．

---

変数変換

---

\begin{align} x &= \mathrm{e}^{\xi} \cos\theta\\ y &= \mathrm{e}^{\xi} \sin\theta \end{align}

(円柱に近いほど粗くなる格子を設定している．)

まず普通の極座標変換$ x = r \cos\theta, \: y = r \sin\theta$をやってから，$ r = \mathrm{e}^{\xi}$, $\displaystyle \frac{\partial}{\partial r} = \frac{\partial \xi}{\partial r} \frac{\partial}{\partial \xi} = \mathrm{e}^{-\xi} \frac{\partial}{\partial \xi}$の変換を行えばいい．

半径方向の速度を$v_r$，動径方向の速度を$v_\theta$とする．すなわち，$ \boldsymbol{v} = v_r \vec{e}_r + v_\theta \vec{e}_\theta$としておく．

(1)式

\begin{align} \nabla \cdot \vec{v} &= \frac{1}{r} \frac{\partial}{\partial r} (r v_r) + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} \\ &= \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} \\ &= \mathrm{e}^{-\xi} \left( \frac{\partial v_r}{\partial \xi} + v_r +\frac{\partial v_\theta}{\partial \theta}\right) = 0 \end{align}

(2)式

・$p$の勾配

\begin{align} \nabla p = \frac{\partial p}{\partial r} \vec{e}_r + \frac{1}{r}\frac{\partial p}{\partial \theta}\vec{e}_\theta = \mathrm{e}^{-\xi} \left( \frac{\partial p}{\partial \xi} \vec{e}_r + \frac{\partial p}{\partial \theta}\vec{e}_\theta \right) \end{align}

・$\Delta \vec{v}$

\begin{align} \Delta\vec{v} &= \nabla(\nabla\cdot \vec{v}) - \nabla\times(\nabla\times\vec{v}) \\ &= \mathrm{grad} \left[ \frac{1}{r} \frac{\partial}{\partial r} (r v_r) + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} \right] - \mathrm{rot} \left[ \frac{1}{r}\left\{ \frac{\partial}{\partial r} (r v_\theta) - \frac{\partial v_r}{\partial \theta} \right\} \vec{e}_z \right]\\ &= \frac{\partial}{\partial r}\left[ \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} \right] \vec{e}_r + \frac{1}{r}\frac{\partial}{\partial \theta}\left[ \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} \right]\vec{e}_\theta \\ &\quad -\vec{e}_r \frac{1}{r}\frac{\partial}{\partial \theta}\left[ \frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r} - \frac{1}{r} \frac{\partial v_r}{\partial \theta} \right] + \vec{e}_\theta \frac{\partial}{\partial r}\left[ \frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r} - \frac{1}{r} \frac{\partial v_r}{\partial \theta} \right] \\ &= \vec{e}_r \left( \frac{\partial^2 v_r}{\partial r^2} + \frac{1}{r}\frac{\partial v_r}{\partial r} - \frac{v_r}{r^2} - \frac{1}{r^2}\frac{\partial v_\theta}{\partial \theta} + \frac{1}{r} \frac{\partial^2 v_\theta}{\partial r \partial \theta} - \frac{1}{r}\frac{\partial^2 v_\theta}{\partial r \partial\theta } - \frac{1}{r^2}\frac{\partial v_\theta}{\partial \theta} + \frac{1}{r^2}\frac{\partial^2 v_r}{\partial \theta^2} \right) \\ &\quad+ \vec{e}_\theta \left( \frac{1}{r}\frac{\partial^2 v_r}{\partial r \partial\theta} + \frac{1}{r^2}\frac{\partial v_r}{\partial \theta} + \frac{1}{r^2}\frac{\partial^2 v_\theta}{\partial \theta^2}+ \frac{\partial^2 v_\theta}{\partial r^2} + \frac{1}{r}\frac{\partial v_\theta}{\partial r} - \frac{v_\theta}{r^2} + \frac{1}{r^2}\frac{\partial v_r}{\partial \theta} - \frac{1}{r}\frac{\partial^2 v_r}{\partial r \partial\theta} \right) \\ &= \vec{e}_r \left( \frac{\partial^2 v_r}{\partial r^2} + \frac{1}{r}\frac{\partial v_r}{\partial r} + \frac{1}{r^2}\frac{\partial v_r}{\partial \theta^2} - \frac{v_r}{r^2} - \frac{2}{r^2}\frac{\partial v_\theta}{\partial \theta} \right) \\ &\quad+ \vec{e}_\theta \left( \frac{\partial^2 v_\theta}{\partial r^2} + \frac{1}{r}\frac{\partial v_\theta}{\partial r} + \frac{1}{r^2}\frac{\partial^2 v_\theta}{\partial \theta^2} - \frac{v_\theta}{r^2} + \frac{2}{r^2}\frac{\partial v_r}{\partial \theta} \right) \\ &= \mathrm{e}^{-2\xi} \left( \frac{\partial^2 v_r}{\partial \xi^2} + \frac{\partial^2 v_r}{\partial \theta^2} - v_r- 2\frac{\partial v_\theta}{\partial \theta} \right) \vec{e}_r \\ &\quad+ \mathrm{e}^{-2\xi}\left( \frac{\partial^2 v_\theta}{\partial \xi^2} + \frac{\partial^2 v_\theta}{\partial \theta^2} - v_\theta+ 2\frac{\partial v_r}{\partial \theta} \right)\vec{e}_\theta \end{align}

・移流項

\begin{align} (\vec{v}\cdot\nabla)\vec{v} &= \left( v_r\frac{\partial}{\partial r} + \frac{v_\theta}{r} \frac{\partial}{\partial \theta} \right) (v_r \vec{e}_r + v_\theta \vec{e}_\theta) \end{align}

と書けるが，$\vec{e}_r,\vec{e}_\theta$は$\theta$によって変化するので，この単位ベクトルの動径方向微分を考えなければならない．

\begin{align} \vec{e}_r &= \cos\theta\vec{e}_x + \sin\theta\vec{e}_y \\ \vec{e}_\theta &= -\sin\theta\vec{e}_x + \cos\theta\vec{e}_y \end{align}

より，

\begin{align} \frac{\partial}{\partial \theta} \vec{e}_r &= -\sin\theta\vec{e}_x + \cos\theta\vec{e}_y = \vec{e}_\theta \\ \frac{\partial}{\partial \theta} \vec{e}_\theta &= -\cos\theta\vec{e}_x + \sin\theta\vec{e}_y = -\vec{e}_r \end{align}

だから，これより

\begin{align} (\vec{v}\cdot\nabla)\vec{v} &= \left( v_r\frac{\partial}{\partial r} + \frac{v_\theta}{r} \frac{\partial}{\partial \theta} \right) (v_r \vec{e}_r + v_\theta \vec{e}_\theta) \\ &=v_r\frac{\partial v_r}{\partial r}\vec{e}_r + \frac{v_\theta}{r} \frac{\partial v_r }{\partial \theta}\vec{e}_r + \frac{v_r v_\theta}{r} \frac{\partial \vec{e}_r }{\partial \theta} \\ &\quad + v_r\frac{\partial v_\theta}{\partial r}\vec{e}_\theta + \frac{v_\theta}{r} \frac{\partial v_\theta }{\partial \theta}\vec{e}_\theta + \frac{v_\theta^2}{r} \frac{\partial \vec{e}_\theta }{\partial \theta} \\ &= \vec{e}_r \left( v_r\frac{\partial v_r}{\partial r} + \frac{v_\theta}{r} \frac{\partial v_r }{\partial \theta} - \frac{v_\theta^2}{r} \right) + \vec{e}_\theta \left( v_r\frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r} \frac{\partial v_\theta }{\partial \theta} + \frac{v_r v_\theta}{r} \right) \\ &= \mathrm{e}^{-\xi}\left( v_r\frac{\partial v_r}{\partial \xi} + v_\theta \frac{\partial v_r }{\partial \theta} - v_\theta^2 \right)\vec{e}_r + \mathrm{e}^{-\xi} \left( v_r\frac{\partial v_\theta}{\partial \xi} + v_\theta\frac{\partial v_\theta }{\partial \theta} + v_r v_\theta \right)\vec{e}_\theta \end{align}

と書ける．(移流項の計算は解析力学のEuler-Lanrange方程式を使って求めることもできる．) ということで，これでやっと(2)式が変換できる．

\begin{align} \frac{\partial v_r}{\partial t} + \mathrm{e}^{-\xi}\left( v_r\frac{\partial v_r}{\partial \xi} + v_\theta \frac{\partial v_r }{\partial \theta} - v_\theta^2 \right) &= -\mathrm{e}^{-\xi} \frac{\partial p}{\partial \xi} + \frac{\mathrm{e}^{-2\xi}}{Re} \left( \frac{\partial^2 v_r}{\partial \xi^2} + \frac{\partial^2 v_r}{\partial \theta^2} - v_r- 2\frac{\partial v_\theta}{\partial \theta} \right) \\ \frac{\partial v_\theta}{\partial t} + \mathrm{e}^{-\xi} \left( v_r\frac{\partial v_\theta}{\partial \xi} + v_\theta\frac{\partial v_\theta }{\partial \theta} + v_r v_\theta \right) &= -\mathrm{e}^{-\xi} \frac{\partial p}{\partial \theta} + \frac{\mathrm{e}^{-2\xi}}{Re} \left( \frac{\partial^2 v_\theta}{\partial \xi^2} + \frac{\partial^2 v_\theta}{\partial \theta^2} - v_\theta+ 2\frac{\partial v_r}{\partial \theta} \right) \end{align}

(3)式

・$ \Delta p$

\begin{equation} \Delta p = \frac{1}{r} \frac{\partial}{\partial r}\left( r\frac{\partial p}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 p}{\partial \theta^2} = \mathrm{e}^{-2\xi} \left( \frac{\partial^2 p}{\partial \xi^2} + \frac{\partial^2 p}{\partial \theta^2} \right) \end{equation}

・移流項の発散

\begin{align} \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) &= \frac{1}{r} \frac{\partial}{\partial r} \left(r [\mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v})_r \right) + \frac{1}{r} \frac{\partial}{\partial \theta} [\mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v})]_\theta \\ &= \frac{1}{r} \frac{\partial}{\partial r} \left(r v_r\frac{\partial v_r}{\partial r} + v_\theta \frac{\partial v_r }{\partial \theta} - v_\theta^2 \right) + \frac{1}{r} \frac{\partial}{\partial \theta} \left( v_r\frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r} \frac{\partial v_\theta }{\partial \theta} + \frac{v_r v_\theta}{r} \right) \\ &= \left( \frac{\partial v_r}{\partial r} \right)^2 + \frac{2}{r} \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial r} + \frac{1}{r^2} \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 \\ & \quad + \frac{v_r}{r} \frac{\partial v_r}{\partial r} + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r} \frac{\partial^2 v_\theta}{\partial r \partial \theta} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \tag{$\star$} \\ &\quad + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r^2} \frac{\partial v_\theta}{\partial \theta^2} + \frac{v_\theta}{r^2} \frac{\partial v_r}{\partial \theta} \tag{$\star\star$} \end{align}

連続の式：

\begin{equation} \nabla\cdot \vec{v} = 0 \Leftrightarrow \frac{1}{r} \frac{\partial}{\partial r} (r v_r) + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} = \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} = 0 \end{equation}

を用いて $(\star), (\star\star)$の部分を変換する．

\begin{align} (\star) &= \frac{v_r}{r} \frac{\partial}{\partial r} \left( v_r + \frac{\partial v_\theta}{\partial \theta} \right) + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ &= \frac{v_r}{r} \frac{\partial}{\partial r} \left( -r \frac{\partial v_r}{\partial r} \right) + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ &= -v_r \frac{\partial^2 v_r}{\partial r^2} - \frac{v_r}{r} \frac{\partial v_r}{\partial r} + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ &= -\frac{v_r}{r} \frac{\partial v_r}{\partial r} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ (\star\star) &= \frac{v_\theta}{r^2} \frac{\partial}{\partial \theta} \left( v_r + \frac{\partial v_\theta}{\partial \theta} \right) + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= \frac{v_\theta}{r^2} \frac{\partial}{\partial \theta} \left( -r \frac{\partial v_r}{\partial r} \right) + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= -\frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial\theta} + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \end{align}

より，

\begin{align} (\star) + (\star\star) &= -\frac{v_r}{r} \frac{\partial v_r}{\partial r} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= \frac{v_r}{r} \left( \frac{v_r}{r} + \frac{1}{r}\frac{\partial v_\theta}{\partial \theta} \right) + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= \frac{v_r^2}{r^2} + \frac{2 v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \end{align}

が得られる．したがって，

\begin{align} \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) &= \left( \frac{\partial v_r}{\partial r} \right)^2 + \frac{2}{r} \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial r} + \frac{1}{r^2} \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 \\ &\quad + \frac{v_r^2}{r^2} + \frac{2 v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ & = \frac{1}{\mathrm{e}^{2\xi}} \left\{ \left( \frac{\partial v_r}{\partial \xi} \right)^2 + 2 \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial \xi} + \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 + v_r^2 + 2\left( v_r \frac{\partial v_\theta}{\partial \theta} - v_\theta \frac{\partial v_\theta}{\partial \xi} \right) \right\} \end{align}

まで変形できた．以上より(3)式は

\begin{align} \frac{\partial^2 p}{\partial \xi^2} + \frac{\partial^2 p}{\partial \theta^2} = -\left\{ \left( \frac{\partial v_r}{\partial \xi} \right)^2 + 2 \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial \xi} + \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 + v_r^2 + 2\left( v_r \frac{\partial v_\theta}{\partial \theta} - v_\theta \frac{\partial v_\theta}{\partial \xi} \right) \right\} \end{align}

というポアソン方程式に変換される。あとは数値計算のために$D$の時間微分の項を付け加えるだけ．

\begin{align} D = \nabla\cdot\vec{v} = \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} = \mathrm{e}^{-\xi} \left( \frac{\partial v_r}{\partial \xi} + v_r + \frac{\partial v_\theta}{\partial \theta} \right) \end{align}

なので，これに$\mathrm{e}^{2\xi}$をかけて$\Delta t$で割ったものを加えることで，数値計算に使う式は，

\begin{align} \frac{\partial^2 p}{\partial \xi^2} + \frac{\partial^2 p}{\partial \theta^2} = &-\left\{ \left( \frac{\partial v_r}{\partial \xi} \right)^2 + 2 \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial \xi} + \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 + v_r^2 + 2\left( v_r \frac{\partial v_\theta}{\partial \theta} - v_\theta \frac{\partial v_\theta}{\partial \xi} \right) \right\}\\ & + \frac{\mathrm{e}^{\xi} }{\Delta t} \left( \frac{\partial v_r}{\partial \xi} + v_r + \frac{\partial v_\theta}{\partial \theta} \right) \end{align}

となる．

渦度

\begin{align} \zeta &= \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \\ &= \left(\cos\theta\frac{\partial v}{\partial r} - \frac{\sin\theta}{r}\frac{\partial v}{\partial \theta}\right) -\left( \sin\theta\frac{\partial u}{\partial r} + \frac{\cos\theta}{r}\frac{\partial u}{\partial \theta} \right)\\ &= \left(\cos\theta\frac{\partial (v_r \sin\theta+v_\theta\cos\theta)}{\partial r} - \frac{\sin\theta}{r}\frac{\partial (v_r \sin\theta+v_\theta\cos\theta)}{\partial \theta}\right) -\left( \sin\theta\frac{\partial (v_r \cos\theta-v_\theta\sin\theta)}{\partial r} + \frac{\cos\theta}{r}\frac{\partial (v_r \cos\theta-v_\theta\sin\theta)}{\partial \theta} \right)\\ &=\frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r} - \frac{1}{r}\frac{\partial v_r}{\partial \theta}\\ &=\frac{1}{e^\xi}\left( \frac{\partial v_\theta}{\partial \xi} + v_\theta - \frac{\partial v_r}{\partial \theta} \right) \end{align}

---

境界条件

---

・$\theta$方向

周期境界

・$\xi$方向

・$\xi=0$：

\begin{equation} v_r=v_\theta=0, \quad \frac{\partial p}{\partial \xi}=\frac{e^{-\xi}}{Re}\frac{\partial^2 v_r}{\partial \xi^2} \end{equation}

・遠方：

$u=1,v=0$すなわち$v_r=\cos\theta, v_\theta=-\sin\theta$、圧力は定数$p=p_0$

---

初期条件

---

\begin{align} v_r&=\left(1-\frac{1}{r^2}\right)\cos\theta \\ v_\theta&=-\left(1+\frac{1}{r^2}\right)\sin\theta \\ p&=p_0-\frac{1}{2}\rho(v_r^2+v_\theta^2) \end{align}

using Plots
Plots.pyplot()

Plots.PyPlotBackend()

i方向が$\xi$，j方向が$\theta$

function BCforP(p)
    for j=1:JM
        p[IM, j]=0.0
        p[1,j]=(4*p[2,j]-p[3,j])/3.0
    end
end

function BCforV(u,v, dtheta)
    for j=1:JM
        u[1,j]=0.0
        v[1,j]=0.0
        
        u[IM,j]=cos((j-1)*dtheta)
        v[IM,j]=-sin((j-1)*dtheta)
    end
end

BCforV (generic function with 1 method)

function SolvePoisson(p, u,v, dtheta, dxi, dt, rhs)
#     rhs=zeros(IM,JM)
    
    for i=2:IM-1
        for j=2:JM-1
            uxi=(u[i+1,j]-u[i-1,j])/(2*dxi)
            vxi=(v[i+1,j]-v[i-1,j])/(2*dxi)
            utheta=(u[i,j+1]-u[i,j-1])/(2*dtheta)
            vtheta=(v[i,j+1]-v[i,j-1])/(2*dtheta)
            
            rhs[i,j] = uxi^2 + 2*vxi*utheta + vtheta^2 + u[i,j]^2 + 2*(u[i,j]*vtheta-v[i,j]*vxi) - 
                    exp((i-1)*dxi)/dt*(uxi+u[i,j]+vtheta)
        end
    end
    
    for i=2:IM-1
        uxi=(u[i+1,1]-u[i-1,1])/(2*dxi)
        vxi=(v[i+1,1]-v[i-1,1])/(2*dxi)
        utheta=(u[i,2]-u[i,JM-1])/(2*dtheta)
        vtheta=(v[i,2]-v[i,JM-1])/(2*dtheta)
            
        rhs[i,1] = uxi^2 + 2*vxi*utheta + vtheta^2 + u[i,1]^2 + 2*(u[i,1]*vtheta-v[i,1]*vxi) - 
                    exp((i-1)*dxi)/dt*(uxi+u[i,1]+vtheta)
    end
    
    @.(rhs[2:IM-1,JM] = copy(rhs[2:IM-1,1]))
    
    for iteration=1:MAXITP
        res=0.0
        for i=2:IM-1
            dp= (dtheta^2*(p[i+1,1]+p[i-1,1]) + dxi^2*(p[i,2]+p[i,JM-1]) + dtheta^2*dxi^2*rhs[i,1])/2.0/(dxi^2+dtheta^2)
            dp=dp-p[i,1]
            res+=dp^2
            p[i,1]+=OMEGA*dp
        end
        p[2:IM-1,JM]=copy(p[2:IM-1,1])
        for i=2:IM-1
            for j=2:JM-1
                dp= (dtheta^2*(p[i+1,j]+p[i-1,j]) + dxi^2*(p[i,j+1]+p[i,j-1]) + dtheta^2*dxi^2*rhs[i,j])/2.0/(dxi^2+dtheta^2)
                dp=dp-p[i,j]
                res+=dp^2
                p[i,j]+=OMEGA*dp
            end
        end
        BCforP(p)
        res=sqrt(res/(IM*JM))
        if res < ERRORP
            resp=res
            itrp=iteration
            break
        end
    end
end

SolvePoisson (generic function with 1 method)

function VelocityEQ(p,u,v,dtheta,dxi,dt,urhs,vrhs)
#     urhs=zeros(IM,JM)
#     vrhs=zeros(IM,JM)
    
    for i=2:IM-1
        urhs[i,1]= -exp(-(i-1)*dxi)*(p[i+1,1]-p[i-1,1])/(2*dxi) + exp(-2*(i-1)*dxi)/RE * ( (u[i+1,1]-2*u[i,1]+u[i-1,1])/(dxi^2) +
                    (u[i,2]-2*u[i,1]+u[i,JM-1])/(dtheta^2) - u[i,1] - (v[i,2]-v[i,JM-1])/dtheta )
        vrhs[i,1]= -exp(-(i-1)*dxi)*(p[i,2]-p[i,JM-1])/(2*dtheta) + exp(-2*(i-1)*dxi)/RE * 
        ( (v[i+1,1]-2*v[i,1]+v[i-1,1])/(dxi^2) + (v[i,2]-2*v[i,1]+v[i,JM-1])/(dtheta^2) - v[i,1] + (u[i,2]-u[i,JM-1])/dtheta )
        
        urhs[i,JM]=urhs[i,1]
        vrhs[i,JM]=vrhs[i,1]
    end
    for i=2:IM-1
        for j=2:JM-1
            urhs[i,j]= -exp(-(i-1)*dxi)*(p[i+1,j]-p[i-1,j])/(2*dxi) + exp(-2*(i-1)*dxi)/RE * 
            ( (u[i+1,j]-2*u[i,j]+u[i-1,j])/(dxi^2) +
            (u[i,j+1]-2*u[i,j]+u[i,j-1])/(dtheta^2) - u[i,j] - (v[i,j+1]-v[i,j-1])/dtheta )
            
            vrhs[i,j]= -exp(-(i-1)*dxi)*(p[i,j+1]-p[i,j-1])/(2*dtheta) + exp(-2*(i-1)*dxi)/RE * 
            ( (v[i+1,j]-2*v[i,j]+v[i-1,j])/(dxi^2) + 
            (v[i,j+1]-2*v[i,j]+v[i,j-1])/(dtheta^2) - v[i,j] + (u[i,j+1]-u[i,j-1])/dtheta )
        end
    end
    
    for i=2:IM-1
        for j=1:JM
            urhs[i,j] -= exp(-(i-1)*dxi) * (-v[i,j]^2)
            vrhs[i,j] -= exp(-(i-1)*dxi) * (v[i,j]*u[i,j])
        end
    end
    
    #移流項 xi方向
    for i=2:IM-1
        for j=1:JM
            if i==2
                urhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j] * (-u[i+2,j]+8.0*(u[i+1,j]-u[i-1,j])+(2.0*u[i-1,j]-u[i,j]))/(12.0*dxi) + 
                abs(u[i, j]) * (u[i+2,j]-4.0*u[i+1,j]+6.0*u[i,j]-4.0*u[i-1,j]+(2.0*u[i-1,j]-u[i,j])) / (4.0*dxi))
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j] * (-v[i+2,j]+8.0*(v[i+1,j]-v[i-1,j])+(2.0*v[i-1,j]-v[i,j]))/(12.0*dxi) + 
                abs(u[i, j]) * (v[i+2,j]-4.0*v[i+1,j]+6.0*v[i,j]-4.0*v[i-1,j]+(2.0*v[i-1,j]-v[i,j])) / (4.0*dxi))
            elseif i==IM-1
                urhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-2.0*u[i+1,j]+u[i,j]+8*(u[i+1,j]-u[i-1,j])+u[i-2,j])/(12.0*dxi) + 
                abs(u[i, j]) * (2.0*u[i+1,j]-u[i,j]-4.0*u[i+1,j]+6.0*u[i,j]-4.0*u[i-1,j]+u[i-2,j])/(4.0*dxi))
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-2.0*v[i+1,j]+v[i,j]+8*(v[i+1,j]-v[i-1,j])+v[i-2,j])/(12.0*dxi) + 
                abs(u[i, j]) * (2.0*v[i+1,j]-v[i,j]-4.0*v[i+1,j]+6.0*v[i,j]-4.0*v[i-1,j]+v[i-2,j])/(4.0*dxi))
            else
                urhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-u[i+2,j]+8*(u[i+1,j]-u[i-1,j])+u[i-2,j])/(12.0*dxi) + 
                abs(u[i, j]) * (u[i+2,j]-4.0*u[i+1,j]+6.0*u[i,j]-4.0*u[i-1,j]+u[i-2,j])/(4.0*dxi))
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-v[i+2,j]+8*(v[i+1,j]-v[i-1,j])+v[i-2,j])/(12.0*dxi) + 
                abs(u[i, j]) * (v[i+2,j]-4.0*v[i+1,j]+6.0*v[i,j]-4.0*v[i-1,j]+v[i-2,j])/(4.0*dxi))
            end
        end
    end
    
    #移流項 theta方向
    for i=2:IM-1
        for j=1:JM
            if j==1
                urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,j+2]+8.0*(u[i,j+1]-u[i,JM-1])+u[i,JM-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (u[i,j+2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,JM-1]+u[i,JM-2]) / (4.0*dtheta) )
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,j+2]+8.0*(v[i,j+1]-v[i,JM-1])+v[i,JM-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (v[i,j+2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,JM-1]+v[i,JM-2]) / (4.0*dtheta) )
            elseif j==2
                urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,j+2]+8.0*(u[i,j+1]-u[i,j-1])+u[i,JM-1])/(12.0*dtheta) + 
                abs(v[i,j]) * (u[i,j+2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,JM-1]) / (4.0*dtheta) )
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,j+2]+8.0*(v[i,j+1]-v[i,j-1])+v[i,JM-1])/(12.0*dtheta) + 
                abs(v[i,j]) * (v[i,j+2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,JM-1]) / (4.0*dtheta) )
            elseif j==JM
                urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,3]+8.0*(u[i,2]-u[i,j-1])+u[i,j-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (u[i,3]-4.0*u[i,2]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,j-2]) / (4.0*dtheta) )
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,3]+8.0*(v[i,2]-v[i,j-1])+v[i,j-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (v[i,3]-4.0*v[i,2]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,j-2]) / (4.0*dtheta) )
            elseif j==JM-1
                urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,2]+8.0*(u[i,j+1]-u[i,j-1])+u[i,j-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (u[i,2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,j-2]) / (4.0*dtheta) )
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,2]+8.0*(v[i,j+1]-v[i,j-1])+v[i,j-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (v[i,2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,j-2]) / (4.0*dtheta) )
            else
                urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,j+2]+8.0*(u[i,j+1]-u[i,j-1])+u[i,j-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (u[i,j+2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,j-2]) / (4.0*dtheta) )
                vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,j+2]+8.0*(v[i,j+1]-v[i,j-1])+v[i,j-2])/(12.0*dtheta) + 
                abs(v[i,j]) * (v[i,j+2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,j-2]) / (4.0*dtheta) )
            end
        end
    end
    
    #時間発展
    for i=2:IM-1
        for j=1:JM
            u[i, j] += dt * urhs[i, j]
            v[i, j] += dt * vrhs[i, j]
        end
    end
end

VelocityEQ (generic function with 1 method)

const IM=81
const JM=121
const MAXITP=10000
const OMEGA=1.0
const ERRORP=1e-4
const MAXSTEP=10000
const RE=100
#const IO_FREQUENCY=50

100

function main(;step=10000,freq=100)
    CFL=0.2

    dxi=0.05
    dtheta=2*pi/(JM-1)
    dt=CFL*min(dxi,dtheta)
    rho=1.0
    p0=0.0

    u=zeros(IM,JM)
    v=zeros(IM,JM)
    p=zeros(IM,JM)
    
    rhs=zeros(IM,JM)
    urhs=zeros(IM,JM)
    vrhs=zeros(IM,JM)
    
    vx=zeros(IM,JM)
    vy=zeros(IM,JM)
    zeta=zeros(IM,JM)
    Ps=Array{Array{Float64,2}}(undef, step÷freq+1)
    Us=Array{Array{Float64,2}}(undef, step÷freq+1)
    Vs=Array{Array{Float64,2}}(undef, step÷freq+1)
    Zs=Array{Array{Float64,2}}(undef, step÷freq+1)

    #初期条件
    for i=1:IM
        for j=1:JM
            u[i,j]=(1.0 - exp(-2*(i-1)*dxi))*cos((j-1)*dtheta)
            v[i,j]=-(1.0 + exp(-2*(i-1)*dxi))*sin((j-1)*dtheta)
            p[i,j]=p0-rho/2.0*(u[i,j]^2 + v[i,j]^2)
        end
    end
    
    BCforP(p)
    BCforV(u,v,dtheta)
    
    num=1
    Ps[num]=copy(p)
    for j=1:JM
        for i=1:IM
            vx[i,j]=u[i,j]*cos((j-1)*dtheta) - v[i,j]*sin((j-1)*dtheta)
            vy[i,j]=u[i,j]*sin((j-1)*dtheta) + v[i,j]*cos((j-1)*dtheta)
        end
    end
    Us[num]=copy(vx)
    Vs[num]=copy(vy)
    
    for j=1:JM
        for i=1:IM-1
            vxp=u[i,ifelse(j==JM,2,  j+1)]
            vxm=u[i,ifelse(j==1,JM-1,j-1)]
            dvx=(vxp-vxm)/(2*dtheta)
            
            vt =v[i,j]
            if i==1
                dvt=(-3*v[1,j]+4*v[2,j]-v[3,j])/(2*dxi)
            else
                dvt=(v[i+1,j]-v[i-1,j])/(2*dxi)
            end
            zeta[i,j]=exp(-(i-1)*dxi)*(dvt+vt-dvx)
        end
    end
    Zs[num]=copy(zeta)
    
    
    # main loop
    for itr=1:step
        SolvePoisson(p, u,v, dtheta, dxi, dt,rhs)
        BCforP(p)
        VelocityEQ(p,u,v,dtheta,dxi,dt,urhs,vrhs)
        BCforV(u,v,dtheta)
        
        if itr%freq==0
            num+=1
            Ps[num] = copy(p)
            for j=1:JM
                for i=1:IM
                    vx[i,j]=u[i,j]*cos((j-1)*dtheta) - v[i,j]*sin((j-1)*dtheta)
                    vy[i,j]=u[i,j]*sin((j-1)*dtheta) + v[i,j]*cos((j-1)*dtheta)
                end
            end
            Us[num]=copy(vx)
            Vs[num]=copy(vy)
            for j=1:JM
                for i=1:IM-1
                    vxp=u[i,ifelse(j==JM,2,  j+1)]
                    vxm=u[i,ifelse(j==1,JM-1,j-1)]
                    dvx=(vxp-vxm)/(2*dtheta)

                    vt =v[i,j]
                    if i==1
                        dvt=(-3*v[1,j]+4*v[2,j]-v[3,j])/(2*dxi)
                    else
                        dvt=(v[i+1,j]-v[i-1,j])/(2*dxi)
                    end
                    zeta[i,j]=exp(-(i-1)*dxi)*(dvt+vt-dvx)
                end
            end
            Zs[num]=copy(zeta)
            print(itr, " ")
        end
    end
    return Ps,Us,Vs,Zs
end

main (generic function with 1 method)

@time Ps,Us,Vs,Zs=main(step=20000,freq=200);

200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 5800 6000 6200 6400 6600 6800 7000 7200 7400 7600 7800 8000 8200 8400 8600 8800 9000 9200 9400 9600 9800 10000 10200 10400 10600 10800 11000 11200 11400 11600 11800 12000 12200 12400 12600 12800 13000 13200 13400 13600 13800 14000 14200 14400 14600 14800 15000 15200 15400 15600 15800 16000 16200 16400 16600 16800 17000 17200 17400 17600 17800 18000 18200 18400 18600 18800 19000 19200 19400 19600 19800 20000  32.324101 seconds (151.83 k allocations: 129.902 MiB, 0.08% gc time)

xi = range(0.0, 4.0, length=IM)
theta = range(0.0, 2*pi, length=JM)
r=exp.(xi)

xx = zeros(IM,JM); yy=zeros(IM,JM)
for i=1:IM
    for j=1:JM
        xx[i,j]=r[i]*cos(theta[j])
        yy[i,j]=r[i]*sin(theta[j])
    end
end

plev=range(-1.0, 0.75, length=31)
ulev=range(-0.5, 1.5, length=26)
vlev=range(-1, 1, length=21)
N=length(Ps)

anim = @animate for t=1:N
    cfp=contourf(xx, yy, Ps[t],levels=plev,
                xlim=(-5,15),
                ylim=(-5,5),
                clim=(-1,0.75),
                color=:plasma,
                aspect_ratio=1,
                title="Pressure",
                cbar=true,clabels=false
                )
    cfu=contourf(xx, yy, Us[t],levels=ulev,
                xlim=(-5,15),
                ylim=(-5,5),
                clim=(-0.5,1.5),
                color=:turbo,
                aspect_ratio=1,
                title="U",
                cbar=true,clabels=false
                )
    cfv=contourf(xx, yy, Vs[t],levels=vlev,
                xlim=(-5,15),
                ylim=(-5,5),
                clim=(-1,1),
                color=:balance,
                aspect_ratio=1,
                title="V",
                cbar=true,clabels=false
                )
    plot(cfp,cfu,cfv, size=(800,800*1.4),layout=(3,1))
end
filename="karman_enchu.gif"
gif(anim, filename, fps=10)

sys:1: UserWarning: The following kwargs were not used by contour: 'label'
[ Info: Saved animation to C:\Users\****\Desktop\Julia\karman\karman_enchu.gif

zlev=range(-5,5,length=41)
contourf(xx, yy, Zs[101],levels=zlev,
                xlim=(-5,15),
                ylim=(-5,5),
                clim=(-5,5),
                color=:curl,
                aspect_ratio=1,
                title="Curl",
                cbar=true,clabels=false
                )

anim = @animate for t=1:N
    cfz=contourf(xx, yy, Zs[t],levels=zlev,
                xlim=(-5,15),
                ylim=(-5,5),
                clim=(-5,5),
                color=:curl,
                aspect_ratio=1,
                title="Curl",
                cbar=true,clabels=false
                )
    plot(cfz, size=(800,800*1.4/3))
end
filename="curl_enchu.gif"
gif(anim, filename, fps=10)

[ Info: Saved animation to C:\Users\****\Desktop\Julia\karman\curl_enchu.gif